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We present a new calculation of the propagation of protons with energies above 10 eV over distances of 
up to several hundred Mpc. The calculation is based on a Monte Carlo approach using the event generator 
SOPHIA for the simulation of hadronic nucleon-photon interactions and a realistic integration of the particle 
trajectories in a random extragalactic magnetic field. Accounting for the proton scattering in the magnetic field 
affects noticeably the nucleon energy as a function of the distance to their source and allows us to give realistic 
predictions on arrival energy, time delay, and arrival angle distributions and correlations as well as secondary 
particle production spectra. 
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I. INTRODUCTION 

The world statistics of ultra high energy cosmic ray 
(UHECR) events of energy above 10 20 eV has now grown 
to 20 events [JjM]. It is very difficult to accelerate parti- 
cles to such high energies in astrophysical shocks, the pro- 
cess thought to be responsible for the majority of the galactic 
cosmic rays [j^]. This has led to a large number of produc- 
tion models, many of them based on exotic particle physics 
scenarios [Q]. The gyroradii of 10 20 eV protons are signif- 
icantly larger than our own Galaxy and this suggests an ex- 
tragalactic origin [@] for any astrophysical scenario (r g = 
lOOkpc x (E/10 2 ^cV) x (lfiG/B) with E and B being the 
proton energy and the magnetic field strength, respectively). 
The large distances between potential UHECR sources and 
Earth leads to another set of problems first pointed out inde- 
pendently by Greisen and by Zatsepin & Kuzmin, now widely 
known as the GZK effect \fy . UHECR protons interact with 
photons of the microwave background radiation and lose their 
energy relatively rapidly during propagation over distances of 
tens of megaparsecs. This should result in a cutoff in the cos- 
mic ray spectrum at an energy just below 10 20 eV. 

Many different calculations [|7[-pj|], performed using vari- 
ous techniques, of the modification of the cosmic ray spectrum 
due to propagation have been published since the original sug- 
gestion. As a result, the general features of the cosmic ray 
spectrum after propagation are well established. Differences 
between the various approaches are, however, significant and 
the accuracy achieved is not sufficient for the interpretation 
of the existing experimental data, and more accurate calcula- 
tions are needed for the expected significant increase of the 
experimental statistics [|l4 - 1 7 1 . 

Previous calculations can be divided into two classes deal- 
ing mainly with: (a) the energy loss processes [0-|l3[], and 
(b) the deflection and scattering of protons in the extragalac- 
tic magnetic field [ ^9| , p0[p4{ | . The first group of calculations 
shows that small differences in the realization of the proton 
energy loss processes generate observable differences in the 



predicted spectra at Earth. Such calculations, however, can- 
not establish an accurate relation between the distance of a 
potential source and the modification of the proton spectrum 
emitted by this source because the influence of the extragalac- 
tic magnetic field is neglected. Among the calculations of the 
second kind, Refs. [|l8|— p0[] do not consider the proton energy 
losses in a satisfactory way, and Refs. [|2lj-[23|] mostly dis- 
cuss their results in a specific context. Only Achterberg et 
al. [ ^4] , ^5| l gi ye a detailed discussion of the fundamental as- 
pects of UHECR propagation in extragalactic magnetic fields, 
which we are interested in here. 

We present here calculations performed with the photopro- 
duction event generator SOPHIA [p6|], which is proven to re- 
produce well the cross section and final state composition in 
nucleon-photon interactions for energies from the particle pro- 
duction threshold up to hundreds of GeV in the center-of-mass 
system. We also account for all other energy loss processes of 
UHECR nucleons, and calculate the proton deflection in the 
extragalactic magnetic field in three dimensions. 

We restrict ourselves to proton injection energies up to 10 22 
eV, and consider (with few exceptions) proton propagation for 
source distances less than 200 Mpc. The calculations are car- 
ried out using a Monte Carlo technique, and we propagate in- 
dividual protons injected as either a mono-energetic beam, or 
with energies sampled from a fixed source energy spectrum. 
This approach has the advantage of representing fluctuations 
in the proton energy losses very well, thereby giving us a good 
handle on the correlations between energy loss, time of flight 
and angular deviation of the flight direction. As we will show, 
these important UHECR characteristics are deeply intercon- 
nected. For a given source distance, there is a strong corre- 
lation between the amount of energy lost, the time delay, and 
the scattering angle. 

Our calculations are thus mainly relevant to scenarios of 
UHECR acceleration at astrophysical shocks, for which 10 22 
eV is a very generous upper energy limit. With this paper we 
wish to establish limits for the distance of potential UHECR 
proton sources as a function of proton energy and the average 
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strength of the extragalactic magnetic field. We also study the 
angular distribution of UHECR with respect to the source di- 
rection (arrival angle) and the time delays after propagation 
over different distances. In addition, the neutrino fluxes pro- 
duced during the propagation are presented. 

The article is organized as follows. We describe the prop- 
agation method, including the relevant features of the event 
generator SOPHIA, in Section 2. Section 3 gives some in- 
teresting results on the propagation of mono-energetic proton 
beams, and compares our results with other work. Section 4 
analyzes the formation and development of the primary and 
secondary particle spectra for protons injected with a power 
law spectrum. In section 5 we discuss the results, present our 
conclusions, and make suggestions for future work. 



II. COSMIC RAY PROPAGATION 

This section provides a description of our simulation code 
for propagating protons in intergalactic space. We treat en- 
ergy losses due to hadronic and electromagnetic interactions 
of the nucleons with photons of the cosmic microwave back- 
ground radiation as well as the deflection of particles by the 
intergalactic magnetic field. Although we present here only 
results on nucleon propagation in random magnetic fields, our 
approach also allows us to follow the particles in complicated 
magnetic field topologies. Because of the time-consuming de- 
tailed simulation of each nucleon propagation path by Monte 
Carlo, the propagation method described below is not suitable 
for calculations involving large cosmological distances. 



A. Interactions and energy loss processes 

Particles of energy E > 10 18 eV interact with photons of 
the cosmic microwave background radiation giving rise to sec- 
ondary particle production and nucleon energy loss. The most 
important processes are: 

• photoproduction of hadrons, and 

• Bethe-Heitler (BH) production of e + e~ pairs by pro- 
tons. 

We also account for the adiabatic losses due to cosmological 
expansion of the Universe, and for the decay of neutrons pro- 
duced in hadronic production process. Since we restrict our 
calculation to models of UHECR acceleration in astrophysi- 
cal shocks, and energies below 10 22 eV, we consider only in- 
teractions with cosmic microwave background photons. The 
calculation of nucleon propagation at higher energies would 
require the use of models of the radio background (see e.g. 
Ref. [27 1). Since we are not presenting results on the develop- 
ment of electromagnetic cascades initiated by secondary parti- 
cles produced in proton-photon interactions, we can safely ne- 
glect interactions on the universal optical/infrared background 
as well. We keep track, however, of the individual energies of 
all secondaries of photoproduction interactions and are thus 



able to show the spectra of neutrinos generated by primary 
protons after propagation over different distances. 

Hadron production and energy loss in nucleon-photon in- 
teractions is simulated with the event generator SOPHIA y26[]. 
This event generator samples collisions of nucleons with pho- 
tons from isotropic thermal or power law energy distributions, 
using standard Monte Carlo techniques. In this paper the code 
has been used with a blackbody spectrum with T = 2.726 
K to represent the cosmic microwave background. Accord- 
ing to the respective partial cross sections, which have been 
parametrized using all available accelerator data, it invokes 
an interaction either via baryon resonance excitation, one- 
particle t-channel exchange (direct one-particle production), 
diffractive particle production and (non-diffractive) multipar- 
ticle production using string fragmentation. The distribution 
and momenta of the final state particles are calculated from 
their branching ratios and interaction kinematics in the center- 
of-mass frame, and the particle energies and angles in the lab. 
frame are calculated by Lorentz transformations. The decay 
of all unstable particles except for neutrons is treated subse- 
quently using standard Monte Carlo methods of particle decay 
according to the available phase space. The neutron decay 
is implemented separately into the present propagation code. 
The SOPHIA event generator has been tested and shown to be 
in good agreement with available accelerator data. A detailed 
description of the code including the sampling methods, the 
interaction physics used, and the performed tests can be found 
in Ref. 

The Monte Carlo treatment of photoproduction is very im- 
portant, because nucleons lose a large fraction of their energy 
in each interaction. As early as 1985 Hill & Schramm [^[j 
pointed out that the use of a continuous energy loss approxi- 
mation for this process neglects the intrinsic spread of arrival 
energies due to the variation of the energy loss AE per inter- 
action, and the Poissonian distribution in the number of pion 
production interactions during propagation. This results in a 
certain "survival probability" of cosmic rays arriving at Earth 
with energies above the GZK-cutoff, as estimated in the as- 
sumption of continuous energy loss. 

Fig. |l]a shows the energy dependence of all parameters rele- 
vant to the average proton energy loss in the microwave back- 
ground (T=2.726 K) for redshift z = 0. The photoproduction 
interaction length A p h for protons is shown as a dashed line. 
Denoting the proton-photon center-of-mass energy by y/s, the 
interaction length can be written as [|l2| 
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Here E (e) is the proton (photon) energy and the proton and 
neutral pion masses are m p and m T o, respectively. The CMB 
photon density is given by n(e) in units of cm -3 eV _1 and 
the photoproduction cross section, a pl (s), is taken from the 
parametrization implemented in SOPHIA. 

The mean energy loss distance x\ OBS (E), shown in Fig. |l]a 
as triple-dot-dashed curve, is calculated as 



x\ oss (E) = 



E 
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The mean energy loss of the nucleon due to the hadron pro- 
duction, (AE), has been calculated by simulating 10 4 inter- 
actions for each given proton energy, resulting in a statistical 
error of the order of 1%. For E > 10 20 eV losses through 
photomeson production dominate with a loss distance of about 
15 Mpc at E > 8 x 10 20 eV. Below this energy, Bethe-Heitler 
pair production and adiabatic losses due to the cosmological 
expansion in the Hubble flow determine the proton energy 
losses. 

Both the photoproduction interaction and the pair produc- 
tion are characterized by strongly energy dependent cross sec- 
tions and threshold effects. Fig. [l]a shows A p h decreasing 
by more than three orders of magnitude for a proton energy 
increasing by a factor of three. After the minimum A p h is 
reached, the proton energy loss distance is approximately con- 
stant. It is worth noting that the threshold region of A p h is very 
important for the shape of the propagated proton spectrum. As 
pointed out by Berezinsky & Grigoreva rfq], a pile-up of pro- 
tons will be formed at the intersection of the photoproduction 
and pair production energy loss distances. Another, smaller 
pile-up will develop at the intersection of the pair production 
and adiabatic loss functions. 

In the current calculation we treat pair production as a con- 
tinuous loss process which is justified considering its small 
inelasticity of 2m e /m p w 10~ 3 (with m e ,m p being the 
electron and proton masses, respectively) compared to pion- 
photoproduction (k w 0.2 — 0.5). We use the analytical fit 
functions given by Chodorowsky et al. [ ]28| ] to calculate the 
mean energy loss distance for Bethe-Heitler pair production. 
This result is in excellent agreement with results obtained by 
simulating this process via Monte Carlo as done by Protheroe 
& Johnson [Q. 

The turning point from pion production loss dominance to 
pair production loss dominance lies at E « 6 x 10 19 eV, with 
a mean energy loss distance of ss 1 Gpc. The minimum of 
the pair production loss length is reached at E m (2 — 4) x 
10 19 eV. For E < (2 - 3) x 10 18 eV continuous losses due 
to the expansion of the universe dominate. For an Einstein-de 
Sitter (flat, matter-dominated) universe as considered here, the 
cosmological energy loss distance scales with redshift z as 

x lossM {E, z) = -J-(l + z)- 3 / 2 » 4000 Mpc (1 + z)~ 3 / 2 , 

(7) 
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FIG. 1. a) Mean energy loss length due to adiabatic expan- 
sion (upper dotted curve), Bethe-Heitler pair production (dash-dotted 
curve), hadron production (triple-dot-dashed curve). Also shown are 
the hadron interaction length (dashed curve) and the neutron decay 
length (lower dotted curve). The solid line shows the total ii oss . 
b) Ratio of mean energy loss length as calculated in Refs. [Js|] (dot- 
ted), [jujj (long-dashed), [jjfshort-dashed), (dash-dotted), Q 
(dashed-dot-dot-dot), and |J25| (thin solid) to the loss length of the 
present work presented in the upper panel. 



for a Hubble constant of Hq = 75 km/s/Mpc, which we 
use throughout this paper. All other energy loss distances, 
£ioss,BH for Bethe-Heitler pair production and xi OS s,ph for 
photomeson production, scale as 



Xi oss (E,z) = (1 + z)- 3 x loss [(l + z)E, z = 0] 



(8) 



We also show the mean decay distance of ~ 9 x 10~ 9 7„ kpc 
for neutrons, where 7„ is the Lorentz factor of the neutron. 
Obviously, neutrons of energy below 10 21 eV tend to decay, 
whereas at higher energies neutrons tend to interact. 

Since the details of the proton energy loss directly affect the 
proton spectra after propagation, we present the ratio of the 
loss distance in previous calculations to that of our work on a 
linear scale in Fig. |l]b. Generally all values of the energy loss 
distance are in a good qualitative agreement. Rachen & Bier- 
mann [ |To| ] treat both Bethe-Heitler and pion production losses 
very similarly to our work except for the threshold region of 
pion production. In the pair production region our work is 
also in perfect agreement with Protheroe & Johnson [ |l2[ ] . An 
overestimate of the loss distance due to pion production of 
- 10 - 20% in Ref. [Q, however, will result in a small 
shift of the GZK cutoff to higher energies in comparison to the 
present calculations. Berezinsky & Grigoreva [^J used a very 
good approximation for the pion production losses, but under- 
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estimate the energy loss in pair production interactions by at 
least 30-40%. The largest deviation of the combined loss dis- 
tance from our model appears in the calculations of Yoshida 
& Teshima |^]. As already pointed out in Ref. [|l2] the largest 
difference occurs at sw 5 x 10 19 eV where Ref. [|] underes- 
timates pair production losses and uses iei oss values larger by 
about 60%, while photoproduction losses are overestimated 
by up to 50%. In the work of Lee [ |l3| ] pion as well as pair 
production losses are treated in fair agreement with our work, 
with differences up to 40% in the threshold region of pion 
production, and 10-20% otherwise. The energy loss code of 
Lee was also used by Sigl and collaborators [22 23]]. The sim- 
ple analytical estimate of photoproduction losses in the recent 
work of Achterberg et al. j p4] , p5| ] underestimates the photo- 
production loss distance by 10—40%, while x\ oss due to pair 
production losses is overestimated by about 20%. 



B. Method of particle propagation 



UHECR propagation involves two main distance scales: (a) 
the hadronic interaction length A p h of typically 3 to 7 Mpc, 
and (b) the much smaller length scale £ mag of typically 10 
kpc needed for a precise numerical integration of the equa- 
tions of motion in a random magnetic field. A straightforward 
Monte Carlo treatment of the propagation using a step size of 
< mag for both hadronic interactions and the equations of mo- 
tion leads to severe efficiency problems for total propagation 
distances of hundreds of Mpc. Hence, the Monte Carlo simu- 
lation is done in the following way. First the path length X^ lst 
from the current particle position to the next possible hadronic 
interaction is determined from 



X&ist — — Aph.min m (£) j 



(9) 



where A p h, m in is the minimum interaction length for hadronic 
interactions (at maximum redshift possible for a given total 
propagation distance) and £ is a random number uniformly 
distributed in (0, 1]. The nucleon is then propagated over the 
path length X&ist in steps of £ mag , and for charged particles 
Bethe-Heitler losses are taken into account and the deflection 
angle is calculated. A hadronic interaction is then simulated 
with the probability A p h : mm/A p h(-E, z), X p h(E, z) being the 
interaction length for the energy E and redshift z. It is shown 
in Appendix A that this method corresponds exactly to a prop- 
agation simulation using Eq. (Q) with A p h(-E, z) for the calcu- 
lation of the interaction distance at each step with the length 
f 

The reduction of the proton energy due to BH pair produc- 
tion and of all nucleons due to adiabatic expansion is cal- 
culated at every propagation step, whereas the correspond- 
ing loss lengths are updated after a simulated path length of 
A p h,tnin an d every photoproduction interaction. In the case of 
neutrons the decay path length is sampled using Eq. (^) with 
the neutron decay length. The smaller of both the hadronic 
interaction and the decay lengths determines then the larger 
scale of the simulation. 



If a photoproduction interaction has occurred, the new en- 
ergy of the proton (neutron) is substituted for the old one, and 
the energies and particle types of the secondary particles are 
recorded. The event generator SOPHIA generates the full set 
of secondary particles, including nucleon-antinucleon pairs. 
Thus the total flux of nucleons after propagation is slightly 
higher than the injected proton flux. Although this is not es- 
sential for the main results of this paper, it may occasionally 
affect the normalization of the proton arrival spectra. 

The propagation is completed when the distance between 
the injection point and the particle location exceeds the pre- 
defined source distance. To obtain precise results for the time 
delay (e.g. total nucleon path length compared to the path 
length of a light ray), the last integration step is adjusted to 
end exactly at the desired distance. 

Particles are injected at a point in space with a randomly 
chosen small angular deviation from the z-axis which defines 
the main propagation direction. The space along the z-axis is 
subdivided into 32 x 32 x 512 cubes of side 250 kpc, each 
filled with a random magnetic field of average strength (B) 
= 10 -9 Gauss (1 nG) [ |29[ 1 satisfying a Kolmogorov spectrum 
with three logarithmic scales. In practice three field vectors of 
random orientation are sampled at scales £ = 1000, 500, and 
250 kpc with amplitudes proportional to 1 : / 3 (see Appendix 
B). The final magnetic field in each of the 250 kpc cubes is the 
vectorial sum of these three vectors. Cyclic boundary condi- 
tions are imposed in case a particle leaves the space of pre- 
calculated magnetic fields. This means that the magnetic field 
experienced by a particle at location x is the same as the field 
calculated at x', 



NiRi, 



x,y, 



(10) 



with Ri being the size of the pre-calculated magnetic field 
region in direction i. Ni is the largest integer number satisfy- 
ing Xi — NiRi > 0. The magnetic field values are refreshed 
after the calculation of 100 propagations to exclude system- 
atic effects by our choice of field vectors. We have verified 
numerically that the magnetic field constructed in this way 
obeys approximately div(B) = and that recalculations of the 
field at smaller intervals do not change the final result. We 
assume that the magnetic field strength does not scale with 
redshift. More information about the implementation of the 
random magnetic field is given in Appendix B. 

The value chosen for £ mag , in principle, depends strongly 
on the average magnetic field and nucleon energy, and is 
a compromise between the precision of the calculation and 
computing time limits. We have chosen £ mas =10 kpc for 
(B) = 1 nG, with an inverse linear scaling for other B values. 
A step size of 1 kpc has been used for short distance prop- 
agations to ensure accurate results for arrival angle and time 
delay distributions. 

Finally, it should be mentioned that the calculation of the 
redshift at a given distance can be done only approximately. 
The reason is the unknown total travel time of a particle from 
the source to Earth at injection time. The actual travel time 
(path length) can be significantly larger than the light travel 
time along a geodesic and is, in general, different for each sim- 
ulated particle trajectory. In the following we use the proper 
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distance-redshift relation to define the redshift of the source 
and along the travel path at observation time. This approxima- 
tion does not strongly affect our results since we consider here 
mainly distances with redshifts smaller than 0.06 and weak 
magnetic fields. However, it should be noted that, in the case 
of a strong magnetic field, cosmological evolution might be- 
come important already at relatively short distances. 



III. RESULTS AND COMPARISON WITH PREVIOUS 
WORK 



In this section, we present results from the simulation of 
proton propagation. We start with mono-energetic proton 
fluxes for which we can compare our results with previous 
work, and which reflect more directly the different treatments 
of the energy loss processes. We then compare results for the 
propagation of protons injected with a power law spectrum. 

One can divide previous calculations into two general 
groups: Monte Carlo based methods, like our own one, and 
analytical/numerical calculations. Protheroe & Johnson [ |l2| ] 
have used a matrix technique to follow the particles over cos- 
mological distances and calculate the 7-ray, neutrino and nu- 
cleon spectra arriving at Earth. The energy loss matrices for 
all particles are calculated with Monte Carlo event genera- 
tors. We have compared our SOPHIA event generator with 
the one of Ref. [ |T2| ] by propagating with the same method an 
E~ 2 proton spectrum with different exponential cutoffs (see 
Eq. ([TT])). For this purpose we have used SOPHIA and the 
event generator of Ref. [12| to calculate the corresponding 
photoproduction matrices and have applied the two matrices 
to propagation over the same set of distances. A compari- 
son of the resulting secondary particle spectra yields excellent 
agreement, pointing to a similar treatment of the particle pro- 
duction process in the different codes. We have also compared 
the matrix method with our Monte Carlo approach by propa- 
gating an exponentially modified power law injection spec- 
trum over 200 Mpc. Again good agreement is found for the 
resulting f M -spectra, while the v e - and neutron spectra are at 
variance with our calculations, which we attribute to a differ- 
ent treatment of the neutron decay. Also, our Monte Carlo 
method results in more losses due to pair production for dis- 
tances > 200 Mpc and a sharp spike at the injection energy for 
very short distance propagation, a consequence of the Poisson 
nature of photon-proton encounters. This feature is discussed 
in detail in Sect. III.A. 

The approach used by Berezinsky & Grigoreva ^ and 
Rachen & Biermann [ |l0| ] is to solve the transport equation 
quasi-analytically by approximating the collisional terms as 
continuous energy loss terms. This does not take into ac- 
count the Poissonian nature of the pion production process as 
pointed out above, and introduces artifacts into the resulting 
nucleon spectra in form of sharp pile-ups. Lee Jl3[ | used a nu- 
merical technique to solve the transport equation for particle 
propagation without using the continuous loss approximation. 

The common assumption in all this work is to consider the 
spatial propagation as strictly along a null-geodesic, with the 



consequence of not being able to gain knowledge about time 
delays and arrival angles of the cosmic rays with respect to 
light and neutrino propagation. 

A hybrid model, combining a Monte Carlo particle trans- 
port code with analytical techniques was presented by Achter- 
berg et al. [M]. Besides simplifying the properties of the 
energy losses by analytical estimates (see Fig. |ljb), this code 
also describes the scattering in the magnetic field as a diffu- 
sion process employing stochastic differential equations. This 
approach has the advantage to allow large propagation steps, 
and is thus computationally very fast, but has a disadvantage 
at small propagation distances which we discuss further be- 
low. Our approach is to use the Monte Carlo technique for 
simulating particle production and to follow closely cosmic 
ray orbits in 3D-magnetic field configurations while traveling 
through the nearby Universe to Earth. This concept, while be- 
ing the most accurate one, limits our propagation calculation 
to small source distances. 



A. Propagation of mono-energetic protons 

In this section we present distributions of arrival energy, 
arrival direction and time delay of the nucleons, as well as 
neutrino spectra, for mono-energetic injection of protons at 
distances of 2, 8, 32, 128 and 512 Mpc from Earth. Protons 
are injected with energy 10 215 eV. At this energy, propagated 
protons can easily suffer several photoproduction interactions, 
and this tends to emphasize the pion production features. 

Fig. ^ shows the distribution of arrival energy of protons 
and neutrons. Clearly visible is the effect of the statistical na- 
ture of photon-proton encounters, also found qualitatively in 
Ref. [Q. At a distance of 2 Mpc, roughly 60% of all injected 
particles do not interact, and this generates a sharp spike at 
the injection energy. This effect due to Poisson statistics re- 
mains visible for distances up to ~ 30 Mpc, showing up as 
a high-energy spike in the cosmic ray spectrum. At larger 
distances, essentially all injected particles undergo interac- 
tions, and therefore, the high-energy spike vanishes. The ar- 
rival energy distributions then become much narrower, and in 
propagation over larger distances would scale simply with the 
energy loss distance for pair production and adiabatic losses, 
modified by the increasing scattering in the magnetic field. 

Fig. ^ shows the distribution of the average time delay of 
the cosmic rays arriving at Earth with respect to propaga- 
tion along a geodesic with the speed of light. This delay is 
caused by scattering of the charged particles by the intergalac- 
tic magnetic field, leading to an increase of the particle's ef- 
fective path length. Thus, the average time delay increases 
with propagation distance, as visible in Fig. ||. Like the ar- 
rival energy distributions, the distributions of the time delay 
also show signs of Poisson statistics, visible especially when 
propagating over short distances. 

The time delay effectively reflects the arrival energy dis- 
tribution idd oc 1 / E\ xx as a result of the random walk pro- 
cess [ p0| , p4| ]. This also emphasizes the importance of an ac- 
curate treatment of energy losses. For example, a direct com- 
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FIG. 2. Arrival spectrum at Earth for mono-energetic injection 
of protons of energy E = 10 21 ' 5 eV and for various source dis- 
tances as indicated. The sharp spike at injection energy for distances 
D < 32 Mpc is due to the low interaction probability within the 
short distance. 



FIG. 3. Time delay of protons injected at different source dis- 
tances and propagated through a random magnetic field of 1 nG. The 
time delay is defined as the propagation time of a particle minus the 
travel time of a light ray along a geodesic. 



parison with the propagation code of Achterberg et al. p4| , [25| ] 
for (almost) the same propagation parameters has shown dif- 
ferences in the time delay up to one order of magnitude for 
D = 32 Mpc. For the same propagation distance, the code 
by Achterberg et al. produces a peak in the arrival energy dis- 
tribution about a factor of 2 lower than found in the present 
work, due to its 20% overestimation of energy losses in the 
photoproduction regime. Together with a difference in the 
magnetic field sampling, which leads to an effective corre- 
lation length £ COII w 390 kpc for the Kolmogorov spec- 
trum used in the present work (see Appendix B) compared 
to £ colT = 1 Mpc for the homogeneous cell approach used 
in Ref. ^M, the observed differences can then be fully under- 
stood by the relation t^ei oc i co „/E 2 rr , as derived in Ref. p4p. 

Protons with injection energy < 10 19 eV suffer mainly con- 
tinuous BH pair production and adiabatic losses that are pro- 
portional to their path length. The substantial deflection in the 
random magnetic field at such energies results in a significant 
increase of the path length. For protons injected at a suffi- 
ciently large distance this can also lead to excessive time de- 
lays. For example, cosmic rays with energy of about 10 19 eV, 
injected at distances greater than 500 Mpc in a 1 nG magnetic 
field, show a time delay exceeding the Hubble time. This gives 
a strict constraint on the cosmic ray horizon. 

The diffusion coefficient for an effective description of the 
scattering process in the magnetic field is strongly energy de- 
pendent, and so is the time delay, idei- To emphasize this cor- 
relation, and demonstrate the advantages of the Monte Carlo 
approach, we show in Fig. Q the scatter plot of proton energy 
versus delay after propagation over 32 Mpc. There is a strong 
correlation suggesting that energy changes of one and a half 
orders of magnitude lead to differences in delay times of more 



than three orders of magnitude, i.e. we find an energy depen- 
dence similar to (idel) °c (B D/E) 2 as derived by Achterberg 
et al. [24] in the small scattering-angle approximation and the 
quasi-linear approximation of wave-particle interactions. The 
correlation becomes less pronounced when propagating over 
significantly larger distances simply because the arrival en- 
ergy distributions become much narrower and the statistical 
nature of the energy loss is smoothed by the prevailing pair 
production and adiabatic losses. This correlation, however, 
would have very important implications for specific models of 
UHECR production, where the duration of an active phase of 
the source competes with the time delay of the protons during 
propagation. The extreme case would be the acceleration of 
UHECR in gamma ray bursts. The particles with the highest 
energies are expected to arrive first, followed by a dissipat- 
ing widening halo of lower energy protons, as emphasized by 
Waxman & Miralda-Escude rfTJ]. 

For large propagation distances, even protons injected with 
10 21 ' 5 eV show time delays that are a considerable fraction of 
the light propagation time (5 — 10% for 512 Mpc). This would 
lead to a limiting proton horizon for a large set of source dis- 
tances and magnetic field values [pij]. 512 Mpc is already a 
limiting horizon for protons injected with 10 19 eV in 1 nG 
fields, as noted above. 

The scattering that leads to time delay also causes angular 
deviations from the direction to the source, as shown in Fig. [| 
for the injection of mono-energetic protons at the same set of 
distances. Note that in our propagation code the 'observer' 
sits on a sphere surrounding the injection point. The angle 
shown is the angle between the particle's arrival direction and 
direction to the injection point. This 'arrival angle' is some- 
what different from the angle between particle's arrival direc- 
tion and the injection direction. This method may lead to an 
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FIG. 4. Scatter plot of time delay versus energy for protons in- 
jected with energy of 10 21 ' 5 eV after propagation over 32 Mpc in 
random B field of 1 nG. 



underestimate of the scattering angle and the time delay when 
the particle fluxes become nearly isotropic and many particles 
have a high probability to scatter back through the 'observer's 
sphere'. It will not, however, affect strongly the results pre- 
sented in this paper, because, as Fig. || demonstrates, we do 
not reach the limit of isotropic 3D diffusion. 

The features of the angular distribution closely follow the 
time delay distributions already shown. For large propaga- 
tion distances, the cosmic ray arrival directions are distributed 
uniformly up to a maximum deflection angle, which increases 
with propagation distance to reach more than 20° at 5 12 Mpc. 
At propagation distances smaller than ^30 Mpc, thus a few 
times the proton interaction length A p h, a peak at small de- 
flection angles occurs due to the effect of Poisson statistics 
for proton-photon interactions. 

Finally Fig. || shows the electron and muon neutrino spec- 
tra generated by the injection of 10 21 5 eV protons at the same 
set of distances. The muon neutrino spectra develop as a func- 
tion of the proton arrival energy spectra folded with the pho- 
toproduction cross section. The fluxes grow with propagation 
distance, and the maximum neutrino energy shifts to lower 
energy reflecting the decreasing proton energy. The growth 
rate with distance decreases for very large distances, where 
the average proton energy significantly decreases and A p h is 
correspondingly significantly longer. 

Electron neutrino spectra show another, very interesting 
feature, that develops with distance. At a minimum distance 
of 2 Mpc the zy e -flux reaches its maximum of 1/2 of the 
spectrum and shows a somewhat wider energy spectrum, en- 
hanced at low energy. At larger distances an additional v e 
component develops at significantly lower energy. As already 
noted in Ref. [g], these are P e 's from neutron decay. The re- 
sulting protons from the decay process carry most of the en- 
ergy, leaving for the P e 's an average energy of only 5 x 10 -4 
of the original neutron energy, and the ^ e -peak is placed at 
about two orders of magnitude to lower energy with respect 




arrival angle <x arp deg 

FIG. 5. Angular distribution of the arrival angle at Earth for 
mono-energetic injection of protons of energy E = 10 21 ' 5 eV, and 
for various source distances as indicated. The magnetic field is 1 nG. 



to the z^-peak. The strength of this component increases with 
distance relative to the direct v e component from ^ decay. 



B. Cosmological modification of the cosmic ray source spectrum 

Berezinsky & Grigoreva [||| introduced the modification 
factor M (E, z) to represent the cosmological evolution of the 
UHECR spectra. M(E, z) gives the ratio of propagated to 
injected protons at the same energy E, for a fixed injection 
spectrum, as a function of the redshift of the injection dis- 
tance compensating for the proton adiabatic losses. M(E, z) 
is thus exactly unity for proton energies below the p7-particle 
production energy threshold. 

At the highest injection energies the modification factor 
shows the GZK cutoff, followed by a pile-up at the crossover 
of photoproduction and pair production energy loss. This 
pile-up is a direct consequence of the resonance nature of 
photoproduction and the hadronic particle production thresh- 
old. The next feature at still lower energy is a shallow dip 
corresponding to the pair production loss, followed by a small 
pile-up below it. The magnitude of the pile-ups and dips de- 
pend not only on the distance and the mean loss distance at 
the photoproduction/pair production crossover, but also on the 
shape of the proton injection spectrum. Flatter spectra create 
bigger pile-ups, because of the increased number of higher 
energy protons that have interacted to lose energy. The pile- 
up energy is linked to the energy where losses due to pair pro- 
duction take over from pion production losses, and is therefore 
strongly dependent on the details of the loss processes in the 
simulations. Fig. ^a shows M(E, z) for propagation without 
magnetic field for the sole reason of comparison with previ- 
ous work. An E~ 2 proton spectrum with a sharp cutoff at 
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FIG. 6. + and v E + £/ e -spectra at Earth after propagating 
a mono-energetic proton beam of energy 10 21 ' 5 eV at distances of 
2, 8, 32 and 128 Mpc (from bottom to top) in a 1 nG intergalactic 
magnetic field. 



E c = 3 x 10 20 eV is injected, and we propagate over a dis- 
tance of 256 Mpc in our calculation (solid line) compared to 
Refs. (dotted line, D=2AQ Mpc), ^ (dashed-dotted line, 
D=256 Mpc), [§ (dashed line, L>=228 Mpc, E c = 10 20 eV) 
and p"3| ] (dashed-dot-dot-dot line, D=256 Mpc). There is ex- 
cellent agreement at all energies with the work of Protheroe 
& Johnson fll2|]. The sharp photoproduction peak of Rachen 
& Biermann ||10|] is an artifact coming from their continuous 
loss approximation for pion photoproduction. As noted previ- 
ously, Yoshida & Teshima [g] used a loss curve which shows 
a significant deviation from that used in the present paper, and 
hence their corresponding pile-up height is also larger than 
in our work. We agree with the position of the pile-up of 
Lee [Oh. However, due to an overestimate of the loss rate 
at this energy, the magnitude of the pile-up in this paper is 
smaller than in our model. The dip just below the pile-up is in 
reasonable agreement with all other works. 

Fig. ^b illustrates the effect of scattering in the magnetic 
field by comparing the resulting corresponding modification 
factors. The 'no scattering' curve (dashed line, as in Fig. ^a) 
is much higher than the more realistic 'scattering curve' in 
the energy range between 10 18 and 10 19 eV. The reason is 
that particles in this energy range have considerable time de- 
lays and correspondingly much higher total energy loss in 
pair production interactions. Another consequence of the in- 
creased proton travel time due to scattering is the develop- 
ment of a higher pile-up at about 10 18 eV, corresponding to 
the large number of particles moved to lower energies from 
the region of that dip. Note that simulation of 10 18 eV par- 
ticles in a 1 nG field is at the threshold of our direct Monte 
Carlo approach, and the calculation is not carried to lower 
energy where it might show an additional pile-up content. 
Fig. Up thus demonstrates the importance of the proton scat- 
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FIG. 7. Upper panel: modification factors for propagation over 
a distance of 256 Mpc without magnetic field after injection of a 
E~ 2 proton spectrum with a sharp cutoff at E c = 3 x 10 20 eV. 
This calculation (solid line) is compared to Refs. [ JTo| ] (dotted line, 
D= 240 Mpc), |[l|] (dashed-dotted line, D=256 Mpc), [|| (dashed 
line, D=228 Mpc, E c = 10 20 eV) and ^ (dashed-dot-dot-dot line, 
D=256 Mpc). Lower panel: comparison of the modification fac- 
tor for rectilinear propagation (dashed curve, 'no scattering' curve) 
and for propagation in a 1 nG magnetic field (solid line, 'scattering 
curve') including the effect of scattering. 



tering in the extragalactic magnetic fields for the shape of the 
final spectrum on arrival at Earth. 

It is important to note that the curves shown in Fig. ^b are 
calculated for a source with unlimited lifetime. In addition, 
by construction, energy loss due to cosmological evolution 
does not enter the modification factor M(E, z). Imposing 
a constraint on the source lifetime will change the modifica- 
tion factor considerably for low energies because, for a given 
distance, the time delay due to the scattering in the turbulent 
magnetic field might become comparable to or even exceed 
the source lifetime. 



IV. FORMATION OF THE PRIMARY AND SECONDARY 
PARTICLE SPECTRA DURING PROPAGATION 

To study the development of the primary and secondary par- 
ticle spectra we followed the propagation of protons injected 
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with a E 2 power law spectrum with an exponential cutoff at 

10 21.5 i e 



dN 

He 



AE' 1 exp[-i;/(10 2i - 5 eV)] . 



(11) 



We recorded the spectra after propagation over 10 Mpc in- 
tervals up to a source distance of 200 Mpc. The results of 
this calculation are relevant for models of UHECR accelera- 
tion at astrophysical shock fronts, although the cutoff energy 
adopted in this calculation is fairly high. 10,000 protons were 
injected with a power law spectrum (integral spectral index 7 
= 1) in each of 30 energy bins covering energies from 10 19 
to 10 22 eV, i.e. 10 bins per decade of energy. We did not 
simulate the propagation of lower energy particles, which do 
not experience photoproduction interactions, but followed the 
secondaries down to arbitrary low energies. 

Fig. ^ shows the evolution of the particles injected in the 
highest energy bin 10 219 to 10 22 eV. The size of each rect- 
angle is proportional to the fractional energy distribution after 
propagation over 10, 20, etc., Mpc. The rate of energy degra- 
dation is dramatic. After only 10 Mpc the spectrum of protons 
injected in a 0. 1 logarithmic bin have spread over one and a 
half orders of magnitude. The width of the energy distribu- 
tion increases with the propagation distance up to ^30 Mpc 
and then decreases. Qualitatively this behavior is very simi- 
lar to the calculation of Aharonian & Cronin [Oj], although 
the direct comparison is difficult because of the different ap- 
proach to the calculation. The average behavior of all protons 
injected with energy above about 3 x 10 20 eV is similar, al- 
though the magnitude of the spread decreases — particles of 
energy below 10 20 eV suffer much smaller losses. After prop- 
agation over about 100 Mpc the spectrum shown in Fig. || is 
already final - it is concentrated within roughly 1/2 order of 
magnitude around ~ 8 x 10 19 eV. This energy slowly de- 
creases because of pair production and adiabatic losses during 
propagation over larger distances, but without change in the 
shape of the distribution. 

The lower panel of Fig. || shows the fractional energy car- 
ried by different particles after propagation in terms of the 
total energy of the protons injected with energy spectrum de- 
scribed by Eq. 11 The proton curve, which also includes 
neutrons, always dominates. The energy content in pro- 
tons, however, is only about 50% of that injected for dis- 
tances above 120 Mpc. The rest of the injected energy is 
distributed between the electromagnetic component and neu- 
trinos. Note the difference between the photon (and elec- 
tron) components from photoproduction (long dashed line), 
and from pair production (short dashed line). While the pho- 
toproduction component rises very quickly and changes very 
little after 100 Mpc, the pair production component is almost 
proportional to the distance, as most of the injected protons, 
despite the high threshold of 10 19 eV, have similar pair pro- 
duction losses. At distances of 100 (200) Mpc 51% (43%) 
of the injected power is carried by nucleons, 31% (37%) by 
the electromagnetic component and 18% (20%) by neutrinos. 
The neutrino fluxes will remain at the same level during prop- 
agation over larger distances, and the established energy bal- 
ance will only slightly change as nucleons yield some of their 
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FIG. 8. a) Arrival energy distribution for protons injected with 
energy between 10 21 ' 9 and 10 22 eV after propagation on 10, 20, ... 
200 Mpc. b) Fractional energy contained in nucleons (solid line), 
7-rays from photoproduction (long dashes) and BH pair produc- 
tion (short dashes) for protons injected with the energy spectrum of 
Eq. IllL The dash-dot lines show the fractional energy in muon (long) 
and electron (short) neutrinos and antineutrinos. 



power to the electromagnetic component through pair produc- 
tion. Adiabatic losses will, of course, affect all components in 
the same way. 

In addition to distributing a fraction of the energy of the in- 
jected protons to secondary particles, the propagation changes 
the energy spectrum of protons. The most energetic nucleons 
lose energy fast and are downgraded after a short propagation 
distance. The number of nucleon with energy above 10 21 eV 
decreases by 10%, 50% and 90% from the injected number of 
protons after only 1, 6, and 20 Mpc. The corresponding dis- 
tances for nucleons of energy above 10 20 eV are 10, 40 and 85 
Mpc. The magnitude of these changes emphasizes the impor- 
tance of detection of very high energy particles: for particles 
of energy above 3 x 10 20 eV (same as the highest energy event 
detected by the Fly's Eye [f30|]) these distances are 1, 10 and 30 
Mpc. The rapid absorption of the highest energy cosmic rays 
implies that the horizon of the highest energy protons is very 
small, and increases the energetics requirements for potential 
UHECR sources. 
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V. DISCUSSION, CONCLUSIONS AND OUTLOOK 



The Monte Carlo propagation of ultra high energy protons 
in a random extragalactic magnetic field has obvious advan- 
tages over other approaches to calculations of proton propa- 
gation in the cosmologically nearby Universe. To start with, 
this approach takes fully into account fluctuations in the po- 
sitions of proton interactions, and thus also in the proton en- 
ergy losses and production of secondary particle fluxes. It also 
naturally generates the correlations between the proton's ar- 
rival energy, its time delay, and its angular deviation from the 
source direction. We have also shown that mathematical ap- 
proaches which use a diffusion description of magnetic scat- 
tering, although superior in computational speed, can lead to 
significant systematic errors for propagation distances smaller 
then ~ 100 Mpc. 

These features of the calculation become extremely valu- 
able when applied to specific models of UHECR acceleration, 
especially models that involve a relatively short (compared to 
light travel time and proton time delay) active phase of the 
source. An extreme example for such a model is the GRB 
model for UHECR acceleration. However, other models in- 
volving interacting galaxies or radio galaxies of specific mor- 
phology could also be affected, especially if embedded in re- 
gions of high (random) magnetic field. 

At energies that allow protons to photoproduce, namely 
above 10 20 eV, the energy degradation is extremely rapid. 
This is not very surprising because of the very short photo- 
production interaction length at energies corresponding to the 
maximum cross section - i.e. A p h below 4 Mpc for energies 
between 4xl0 20 eV and 10 21 eV. This energy range is very 
relevant, as it is just above the highest ener gy p articles de- 
tected by the Fly's Eye and AGASA arrays j3C|,p|]. A large 
part of this rapid energy dissipation in our calculation is due 
to the correct implementation of the fluctuations in photopro- 
duction interactions in SOPHIA. A good example for the size 
of the fluctuations is the proton energy distribution after prop- 
agation over 10 Mpc shown in Fig. [| which covers more than 
one and a half orders of magnitude. This is an extreme case. 
However, every particle injected with an energy well above 
the photoproduction threshold would very rapidly result in a 
distribution extending down to the threshold, within the first 
10 Mpc. 

This rapid energy dissipation creates additional problems 
for models of cosmic ray acceleration at astrophysical shocks. 
Apart from the difficult question of the maximum acceleration 
energy, such models require that a significant fraction (0.01 
to 0.1) of their source luminosity contributes to the UHECR 
flux. The rapid energy dissipation increases the energy re- 
quirements in terms of total luminosity and severely limits the 
source distance. Because of magnetic scattering, such limits 
could also be set for particles injected with energy below the 
photoproduction threshold. 

Fig. ^| shows the 50% horizon for UHECR sources as a 
function of source particle energy for (B) values of 0.1, 1 
and 10 nG. The 50% horizon i?5o is defined here as the light 
propagation distance to the source at which 1/e of all injected 
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FIG. 9. Proton 50% horizon as a function of injection energy for 
average random magnetic fields of 0. 1 (dashed histogram), 1 (solid 
histogram), and 10 (dotted histogram) nG. See text for definition. 
The solid line is the total energy loss length from Fig. [jj shown here 
for comparison. 



protons have retained 50% or more of their energy, i.e. R$q is 
achieved when 



dN 
dE 



dE = iVoexp(-l), 



(12) 



where No is the number of particles injected with energy Eq- 
To start with, R 50 is small at any energy, and demonstrates 
the resonant nature of the photoproduction cross section. At 
E = 10 20 eV R 50 is about 100 Mpc, while at 2xl0 20 eV it 
decreases to 20 Mpc and becomes smaller than 10 Mpc for 
energies above 3x 10 20 eV. For injection energies above 10 20 
eV the horizon energy dependence is similar to that of the 
energy loss distance shown in Fig. |l[ These protons are not 
affected much by the magnetic field since their scattering an- 
gles are small, but suffer mainly from energy degradation due 
to %r) encounters. Below 10 20 eV the picture changes. The 
scattering in the magnetic field increases the propagation time 
and thus causes additional energy loss and an increase of the 
ratio x\ oss /Rr M . 

Stronger magnetic fields create delays, that could be longer 
than the light propagation time from the source and reverse 
the trend - the horizon starts decreasing below ~6x 10 19 eV 
and is restricted to 75 Mpc at 10 19 eV. Since the average time 
delay is inversely proportional to E 2 , the decrease of R50 is 
expected to become more drastic at lower energy. One con- 
sequence of the strong energy dependence of R50 is, for ex- 
ample, that our attempts to correlate the arrival directions of 
UHECR with different types of astrophysical objects should 
use only objects within the particle horizon depending on the 
magnetic fields strength in different regions of the Universe. 
Independently of the magnetic field value, however, the hori- 
zon defined above is much smaller than the conventional num- 
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bers of 50 or 100 Mpc for the highest energy cosmic ray 
events. 

There are many relevant astrophysical problems which can 
be studied with the approach described in this paper. We plan 
to use the code for proton propagation in regular magnetic 
fields associated with large scale structures (local superclus- 
ter, supergalactic plane). The regular fields, especially if they 
reach the observationally allowed limits of 0.03 fiG and even 
0.1 //G, could change the propagation patterns for 10 19 eV 
cosmic ray protons and alter the horizon values shown in 
Fig. H We also plan to set limits on models of slow UHECR 
acceleration on shocks of very large dimensions and to look 
for possibilities of ultra-high energy 7-ray halos around the 
sources and along the tracks of the UHECR protons. 
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APPENDIX A: MONTE CARLO SAMPLING OF 
INTERACTION POINTS 

In the following we discuss the application of the veto al- 
gorithm to the sampling of interaction points along a nucleon 
propagation path. The probability of having no hadronic inter- 
action with a photon of the CMB within a path length interval 
(si,s 2 ) reads 

p ^ S2) ^ exp {-JIx^m}- (A1) 

The interaction length itself depends only on the nucleon en- 
ergy. However, because of the treatment of Bethe-Heitler 
losses as continuous process, this energy depends on the path 
length s. Correspondingly, the probability for one interaction 
in the interval (s, s + ds) is given by 

ds 

P int (s)ds = P no (0, s) - (A2) 
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where P no (0, s) is the probability that no interaction has oc- 
curred before. In our approach we replace X p h(E(s)) by the 
constant A p h, m in and use (A.2) to sample the path length dis- 
tance from the current location (s = 0) to the next interac- 
tion. This interaction point is then accepted with the probabil- 
ity A p h,min/Aph(£'(s)). Hence the interaction probability can 
be written as 



P lnt (s)ds 



dsi 



Aph )m i n 



P„o(0,Sl) 1 



A 



ph,min 



ds\ 



^ph,min 



-Pno(0,Sl) 1 



A ph (£(si)) 

^ph,min 

Aph(^(si)) 



dsn 



s± Aph,min 



-Pno(si, S 2 ) 1 - 



A 



ph.min 



+ 



A ph (£(s 2 )) 

Aph,min \ ds 



X ph (E(s))J A ph 



where we have used 



-Pno(S2,Sl) = exp 



A 



ph,min 



-Pno(si,s) 

Pno(s2, s) 
(A3) 

(A4) 



The first term in square brackets corresponds to the probabil- 
ity that no interaction was sampled in the interval (0, s). The 
second term is the contribution which comes from an inter- 
action point sampled at si but rejected with the probability 

1 Aph,min/Aph- 



The integration limits in (A3) ensure the ordering of the 
interaction points according to the simulation method, < 
si < S2 < ■•• < s. Symmetrizing the integration limits 
yields 



P int (s)ds 



ds 



A P h(£(s)) 



exp 
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71 = 



ds' 



^ph,min , 

l l 



exp 



Aph,min 

ds' 1 ds 



o X ph (E(s'))j X ph (E(s)y 



(A5) 



which is identical to ( A2 ) and shows that the described simu- 
lation method reproduces the correct, energy-dependent inter- 
action length. 



APPENDIX B: IMPLEMENTATION OF THE MAGNETIC 
FIELD 

A turbulent magnetic field which is frozen into a fluid with 
fully developed hydrodynamic turbulence would follow a Kol- 
mogorov spectrum, which is defined by 



where fc is the wavenumber [[[l]]. I{k) is the energy density 
per unit wave number, fco the smallest wavenumber of the tur- 
bulence, the inverse k^ 1 is sometimes called the "cell size" 
of the turbulence. Hence we have for the total energy density 
PI 



8tt 



dkl(k) 



(B2) 



In the propagation program we consider 3 discrete wave num- 
bers. Thus we have to rewrite this integral in terms of a 
discrete spectrum in fc, starting with fco an d continuing with 
fej = 2fej_i, i = 1,2. These are equally spaced apart in log 2 k, 
with A(log 2 k)=l. Hence the energy density we should as- 
cribe to each of the three wavenumbers is approximately 



Ui 



I{k) dk 



d(log 2 fc) 
Jo fco m 2 



A(log 2 fc) 

2/3 



fc; 
fc 

The total energy density is then a simple sum, 



8tt 



U + Ui + u 2 



(B3) 



(B4) 



We normalize the field to a total energy density corresponding 
to{\B\) = InG, i.e. U tot « 4x 10" 20 erg cm" 3 . 

The technical implementation of the magnetic field into our 
propagation code is as follows. We divide the propagation 
volume into cubes of 1 Mpc side length, and attach to each 
of them a homogeneous field Bo with magnitude Bq and ran- 
dom direction. Each of these cubes is divided into 8 cubes 
of 0.5 Mpc side length, to which a field Bi of magnitude E>\ 
and random direction is vectorially added to the field Bo- The 
procedure is repeated once more, so that our field is eventually 
realized on elementary cubes of 0.25Mpc side length, each of 
which carries a magnetic field Bo + Bi + B2. We check that 
divB ~ by approximating the surface integral with the sum 
of the outward normal component of B over the surface of the 
8x8x128 Mpc 3 volume V. The volume averaged value of 
div(B) is calculated as 



(V-B) = ^X)B x <fe 



(B5) 



The r.m.s. value of (V • B) for 10,000 field realizations is 
(V ■ B) rms = 3.7 x 10" 6 nG/kpc. 

We also calculate the effective correlation length l COTT by 
equating 



(B(x).B(x + 0) = B 2 (x)exp 



The best fit value of £ COII is 390 kpc. 



lei 



(B6) 



I(k)=I (k/k a )- 5 / 3 



(BI) 
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